Interface Width and Bulk Stability: requirements for the simulation of Deeply 

Quenched Liquid-Gas Systems 
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Simulations of liquid-gas systems with extended interfaces are observed to fail to give accurate 
results for two reasons: the interface can get "stuck" on the lattice or a density overshoot develops 
around the interface. In the first case the bulk densities can take a range of values, dependent on 
the initial conditions. In the second case inaccurate bulk densities are found. In this communication 
we derive the minimum interface width required for the accurate simulation of liquid gas systems 
with a diffuse interface. We demonstrate this criterion for lattice Boltzmann simulations of a van 
der Waals gas. When combining this criterion with predictions for the bulk stability we can predict 
the parameter range that leads to stable and accurate simulation results. This allows us to identify 
parameter ranges leading to high density ratios of over 1000. This is despite the fact that lattice 
Boltzmann simulations of liquid-gas systems were believed to be restricted to modest density ratios 
of less than 20 
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Application of lattice Boltzmann methods to the simu- 
lation of liquid-gas systems has been one of the early suc- 
cessful applications of lattice Boltzmann. Two very dif- 
ferent algorithms were developed to do this: Swift et. al. 
@)Hj0,@ developed an algorithm based on implementing 
a pressure tensor and Shan et al. @, 0, HI developed an 
algorithm based on mimicking microscopic interactions. 
These algorithms have been succesfully applied to sim- 
ulations of phase-separation drop-collisions [3, ITol | . 
wetting dynamics and spreadi ng ll| . and the study of 
dynamic contact angles 12, 12, M | ■ Only recently has it 
been shown that both algorithms perform very simil arly 
when higher order corrections are taken into account [1 51 ] . 

Previously there have been only heuristic analysis of 
what range of paramteres lead to accurate simulation re- 
sults. In some parameter ranges not too far from the 
critical point too thin interefaces can lead to interfaces 
that stick to the lattice. This can lead to non-unique 
bulk densities for the liquid and gas phases that depend 
on the inital contitions. Further from the critical point 
density overshooting is observed at the interfaces. This 
also leads to incorrect bulk densities. In this communica- 
tion we present a criterion that allows us to predict the 
range of acceptable values for the interface width. We 
will show that the accuracy of the algorithm rapidly de- 
teriorates when this limit is exceeded. The second impor- 
tant contribution of this communication is a more general 
definition of the equation of state. Usual lattice Boltz- 
mann methods recover the ideal gas equation of state 
with a pressure of p = p/3 when the gas is dilute. While 
relaxing this requirement has no effect on the interfacial 
properties it allows us to adjust the bulk-stability of the 
lattice Boltzmann method. Taken together this allows 
us to determine sets of paramters for which very deep 
quenches can be simulated. This is an important result 
since lattice Boltzmann methods were previously believed 
to be limited to density ranges of about 20. We can now 
identify parameter ranges for which we can obtain much 



larger density ratios. In this communication we present 
some simulations with density ratios larger than 1000. 

The key to a successful simultion of liquid-gas sys- 
tems is the faithful representation of the interface. In 
particular we need to obtain a constant pressure across 
a flat interface. For a standard Landau Free energy of 
F = / Vo + «/2 (Vp) 2 , where p is the density, we have 
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where po = pd p ipo — ipo is the bulk pressure. In equilib- 
rium the normal pressure is constant across an interface. 
For the equilibrium profile p(x) corresponding to a bulk 
pressure of ps this implies 



Pip) ~ Pb 
pd 2 p - \d x pd x p 
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A problem arises for discrete simulation methods because 
the density derivatives have to be replaced by discrete 
derivatives. For simplicity we will limit our analysis to 
the one-dimensional case here. The standard definitions 
for the discrete derivative are d x p = 0.5[p(x+l)— p(x— 1)] 
and d^p — p(x + 1) + p(x — 1) — 2p(x). This puts severe 
limits on the allowable values for k. 

We can now perform a simple estimate of the minimum 
value K m that allows this pressure to be the equilibrium 
pressure. For any point on the interface with density p s 
we consider two neighboring points, one with a smaller 
density p- and one with a larger density p+. If pi is the 
gas density and p g is the gas density we deduce a lower 
limit for the smallest possible value K m (p s ) by varying 
the values of p+ and p_, but not allowing overshooting: 



K m (p s ) = min — — 
p g <P-<p 3 Ps\P- 
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Since we need to allow any value of p s between p g and pi 
the minimum allowable value of k is then given by 
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To be concrete and to illustrate the power of the min- 
imum k prediction we will consider the pressure of a van 
der Waals gas 
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Please note that we introcuded a scale-factor po with the 
pressure. Such a scale factor clearly is not expected to 
change the equilibrium behavior. However, as we will see 
below, it can have a profound effect on the bulk stability 
of the system. Most previous aproaches correspond to a 
choice of po = 1 in Eq. ((5|) . Only for this choice will the 
lattice Boltzmann method recover the standard lattice 
Boltzmann method for ideal gases in the limit of small 
densities. 

A good approximation for the interface shape that be- 
comes exact close to the critical point is given by 



p mlt {x)= Pa + 
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where pi and p g are the equilibrium gas and liquid den- 
sities. The interface width is given by 
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This profile is not the exact analytical solution to the 
differential equation VP = 0, but it is very close to it. 

We first examine numerically which values of p s in Eq. 
<j3j) lead to the most restrictive constraint, i.e. the largest 
value of K m (p s ). The orange line in the inset of Figure 
[T] shows how this density p cr it varies as a function of 
temperature. Near to the critical temperature, the or- 
ange line in the inset in Figure Q] lies close to the high 
density spinodal curve. This is because P — pb has its 
highest magnitude here, therefore helping to maximize 
K m within this region. An analytical estimate for n m 
can be obtained by expanding the pressure around the 
critical density, giving 
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Within this regime, P—pb is large and negative, therefore 
p_ and p+ must be chosen to make the denominator in 
([3]) as negative as possible. A suitable choice is p_ = p g 
and p+ = p s . We assume that the critical value of p s lies 
on the spinodal curve p S pm = 1 + 2\J6 C — 8. This allows 
us to obtain K m , and substituting this expression into ([7|) 
gives a minimum interface width of 
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Numerical Solution 
a Analytic Solution Low 6/9. 

□ Analytic Solution High 6/9 
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FIG. 1: (Color online) The two limiting cases for which we 
can obtain an analytical approximation to the w(K m ) relation. 
The inset shows the value of the critical density pcrit, which 
is the value p 3 takes when Q is maximized. Note that there 
is a discontinuity. 



As the temperature is decreased in the inset of Figure Q] 
the critical density p cr n makes a discontinuous jump to 
a regime in which it lies close to the gas density, p g . The 
minimum interface width can, in this case, be analytically 
obtained by expanding densities around p g . We define 
p s = p g + dp and p + = p_ + Ap. Since P — pb is a 
positive quantity, a suitable choice for p_ is p_ = p s . 
Substituting these expressions into equation ([3]) gives 
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Minimizing this with respect to Ap leads to p s = Ap/4. 
Re-substituting this result back into equation p0[) . and 
maximizing with respect to Sp, we finally obtain p cr it — 
2p g . Using this we can calculate the minimum interface 
width, 
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as shown by the triangles in Figure[T] This closely follows 
the numerical result at low temperatures. This means 
that we need wide interfaces for deep quenches because 
of the unfortunate cancellation of the discrete derivative 
and Laplace operator for low densities in the denominator 
of©. 

Most previous lattice Boltzmann simulations ap- 
proached the simulation of non-ideal systems by using 
the ideal gas equation of state p = p0 = p/3, as a starting 
point. Interactions are then included to allow the simula- 
tion of non-ideal systems. The speed of sound c s — \Jd p p 



will then recover the ideal gas value of 1/ \/?> in the dilute 
limit. For a van der Waals gas with a critical density of 
1 and a temperature of 9 — 1/3 and an interfacial free 
energy of J re/2 (Vp) 2 the pressure tensor used by previ- 
ous approaches matched the ideal gas equation of state 
in the dilute limit, leading to po = 1. For the van der 
Waals gas the speed of sound increases rapidly for high 
densities. A problem arises when the speed of sound 
becomes larger then the lattice velocity \vi\, because in- 
formation can not be passed on at speeds larger than the 
lattice velocity. When the speed of sound is increased 
above 1 the simulation becomes unstable. This problem 
is exacerbated by the presence of the gradient terms in 
the pressure tensor. These terms further decrease the 
stability, as shown in a previous analysis of the pressure 
method by C. Pooley for one, two, and three dimensional 
lattice Boltzmann methods [17| . In the notation of this 
letter the linear stability condition is 



< V 1 - 4p «P 
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for a homogeneous, one dimensional system with den- 
sity p. This implies a restriction for both the maximum 
quench depth and the maximum intereface witdth. This 
is shown in Fig. [3] as solid lines for different values of po. 
This suggests that, at least as far as the stability of the 
bulk phase is concerned, the most stable solutions should 
be found for k — 0. 

For simplicity we demonstrate the constraints of liquid- 
gas lattice Boltzmann simulations by the common one 
dimensional projection of one, two and three dimensional 
models. This is the D1Q3 model. The lattice Boltzmann 
equation for densites fi corresponding to velocity Vi is 
given by 

/i(x + Vi , t + 1) = /;(x, t) + i(/P(x, t) - /i(x, t)). (13) 

T 

The J® are the equilibrium distribution and are given by 
/_i = -pu/2 + II/2; f = p - II; /i = pu/2 + U/2 



and II = pu 2 + P + vud x p [3, 15 1 where v = (r — 0.5)/3. 
To second order the resulting equations of motion are, as 
usual, the continuity equation 



d t p + d x (pu) = 



(14) 



and the Navier Stokes equation 

d t (pu) + d x {pu 2 ) = -d x P + d x {2vpd xP ) (15) 

To lower the speed of sound in the liquid phase we now 
reduce the value of po in ([5]) . This decreases the speed of 
sound in the liquid by a factor ^Jp~q. This also increases 
the range of stability for re in (fl2"|) . We now expect that 
lowering the speed of sound by a sufficient factor will 
reduce the speed of sound sufficiently to simulate systems 
with arbitrarily low temperature ratios 




FIG. 2: (Color online) The van der Waals phase diagram 
is recovered to very good approximation for interfaces wider 
than the minimum width. The large symbols represent the 
predicted points at which the simulations are predicted to 
become inaccurate by Eq. <[3j) . The value of po only affects 
the bulk stability and values from 1 to 10 -7 were used for 
increasing quench depth. 



To test this idea we performed simulations with near 
equilibrium profiles using a one dimensional three veloc- 
ity Vi = { — 1,0,1} model by defining an initial density 
profile that is given by two domains with densities pi 
and p g respectively connected by an equilibrium interface 
given by Eq. © . This profile is not the exact analytical 
solution to the differential equation VP = 0, but it is 
very close to it. By initializing the simulation with this 
profile we test the linear stability of the method around 
an equilibrium profile to good accuracy. The shape of a 
stable interfacial profile is independent of po- 

In Figure [3] we see that by lowering po the method is 
now able to simulate very small values of the reduced 
temperature 9/9 c for interface width w > 1, but that 
significantly larger widths are required to recover an ac- 
curate phase diagram for deep quenches. For values of 
9/9 c between 0.9 and 1 we also find non unique solutions 
for small values of re which is discussed in more detail in 
a previous paper 
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We now test the predictions of the accuracy of ex- 
tended interface liquid-gas simulations using a lattice 
Botzmann implementation first presented in [151 ]. 

For small values of re the interface becomes sharp in 
the continuous limit so that the derivatives become ar- 
bitrarily large. But in the numerical implementation 
the derivatives are discrete. The discrete values are 
limited by the lattice spacing. In the one dimensional 
case we choose Vp(x) = 0.5(p(x + 1) — p(x — 1)) and 
V 2 p{x) = p(x + 1) - 2p(x) + p(x - 1). 

The methods always lead to a constant pressure, even 
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FIG. 3: (Color online) Existence of accurate solutions for dif- 
ferent values of po and w. Symbols indicate parameter com- 
binations that lead to stable, accurate, and unique solutions. 
Solid lines are the bulk stability limits for the pressure method 
given by Eq. (I12|l . The dashed line is the line for an accurate 
interface representation given by Eq. ([3]). 



across an interface |15j. We performed a scan of the 
parameter space w and 0/9 c initializing the simulation 
with a near equilibrium profile for different values of 
Pq. We accept simulations that are stable, accurate 
and unique. The criterion of accuracy is denned to be 
logio(pmin) - logi (/O s ) < 0.1. As can be seen in Figure 
the results are not very sensitive to the exact value of 
the cutoff. For values of the interface width w < 1.5 we 
also test the uniqueness of the simulation by using initial 
profiles with bulk densities corresponding to the pressure 
at the spinodal points (IBj. Our criterion for uniqueness 
is that all simulations lead to the same minimum density 
to within Ap < 0.01. 

Comparing ([3]) , shown as a dashed line in Figure 02 
and the numerical results for stable, accurate and unique 
solutions shows excellent agreement. The bulk stability 
of Eq. gives the second limit for the acceptable pa- 
rameter range for the pressure method. We perforemd 
a similar analysis for the forcing method of [15| | and ob- 
tained nearly identical results except for a slightly larger 
range of bulk stability. Note that previous lattice Boltz- 
mann simulations use po = 1, which corresponds to the 
area under the black line in Fig. [3] This is why it was 
assumed that lattice Boltzmann simulations are limited 
to a maximum density ratio of about 10 [![. 

The interface constraint ([3]) is remarkably successful 
at predicting the acceptable simulation parameters. It 
predicts how thin is too thin for an interface. It thereby 
detects when non-unique solutions occur and when solu- 
tions for deep quenches fail to deliver accurate results. 

For simulation methods it is important to be aware 
of the acceptable paramter ranges. Lattice Boltzmann 



simulations are often believed to have the nice prop erty 
of becoming unstable before they become inaccurate [18| . 
This is not the case for thin interfaces in multi-phase 
simulations. In this case the simulation can remain stable 
and become inaccurate. This makes it necessary to find 
some criterion that determines whether a set of paramters 
will lead to an accurate simulation. Such a criterion was 
presented for the interface width (controlled by n) in this 
communication. 

The second contribution presented in this paper ap- 
pears trivial at first: in consists of a simple pre- factor for 
the pressure. It is new, because it breaks with the idea 
that the standard lattice Boltzmann method for ideal 
gases should always be recovered in the dilute limit. This 
prefactor has profound implications on the stability of the 
bulk phases, which can be seen using an important result 
about the bulk stability from [l7j |. 

Combining these two components we were able to show 
that lattice Boltzmann is indeed able to simulate very 
deep quenches for liquid-gas cases. This analysis was 
general and will be applied to other equation of states 
as well as other discretization of the interfacial terms. 
This may yet yield significant further advances for the 
development of multi-phase lattice Boltzmann methods. 
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